To determine the composition and structure of an underground formation, seismic data may be used. When a seismic wave propagates through the earth from a source, the seismic wave may propagate until the seismic wave is reflected, such as from an interface. The reflected wave is received by a receiver, the reflected wave may provide information regarding the geophysical properties of the subsurface volume of the underground formation. The properties may include, without limitation, seismic wave velocity, mass density, and anisotropy properties. The time between transmission of the seismic wave and reception of the reflected seismic waves may be measured to determine the distance to an interface. In exploration geophysics, by recording the seismic waves at the surface by one or more geophones or hydrophones, a model of the subsurface strata may be created through various seismic processing techniques.
Computer methods for imaging seismic data may be utilized in oil, gas, and coal exploration to generate maps of the subsurface structure of rock formations. The maps and data may be analyzed to assist in locating oil, gas, or coal, and to identify underground formations.
A conventional technique in computer methods for imaging seismic data is reverse time migration (“RTM”). RTM is a two-way depth migration solution that can accurately describe wave propagation in complex media. However, RTM requires large amounts of computation time and memory.
As a seismic wave propagates through a medium, the elastic energy associated with the wave is gradually absorbed by the medium, resulting in dissipation of the wave as heat energy. Thus, seismic data may be adversely affected by absorption of seismic waves by faults and fracturing in the underground formation and by different solids, fluids, gases and voids. This absorption is known as anelastic attenuation. The effect of anelastic attenuation may be amplitude attenuation (reduction in amplitude) and phase velocity change, i.e., velocity dispersion, including, for instance, frequency dependent phase velocity.
One conventional computer method for correcting for amplitude attenuation and velocity dispersion is “Q-RTM.” The “Q” in “Q-RTM” is the anelastic quality factor. “Q” quantifies the effects of anelastic attenuation. “Q” is inversely proportional to the attenuation. Q may be defined as:
Q=2π(E/δE) where E/δE is the fraction of energy lost per cycle.
In certain embodiments, Q may be constant or nearly constant over a predetermined seismic frequency range.
By introducing memory variables to approximate constant Q effect, a wave propagation equation that incorporates anelastic attenuation may be represented as:
                    ∂        v                    ∂        t              =                  -                  1          ρ                    ⁢              ∇        p                                ∂        p                    ∂        t              =                            -                      λ            ⁡                          (                              1                +                                  L                  ⁢                                                                          ⁢                  τ                                            )                                      ⁢                  ∇                      ·            v                              -                        ∑                      i            =            1                    L                ⁢                  r          i                    +      Q                          ∂                  r          i                            ∂        t              =                            -                      1                          τ                              σ                ⁢                                                                  ⁢                i                                                    ⁢                  r          i                    -                                    λ            ⁢                                                  ⁢            τ                                τ                          σ              ⁢                                                          ⁢              i                                      ⁢                  ∇                      ·            v                              Where t is time, v is the particle velocity, p is the measured pressure of the seismic wave, λ is the bulk modulus, ρ is the mass density, r is the memory variable, S is the source term, τ and τσ are parameters selected to approximate a constant Q over a specific frequency range, and L is the number of relaxation mechanisms.
In certain traditional embodiments, a τ method is designed to simulate anelastic attenuation effects, such as Q attenuation, but cannot be applied to simulate compensation effects, such as Q compensation. From the τ method, an equation for complex velocity has been derived, allowing calculation of Q from the complex velocity equation. Q at frequency ω may be calculated as:
      Q    ⁡          (      ω      )        =            1      +              L        ⁢                                  ⁢        τ            -                        ∑                      l            =            1                    L                ⁢                  τ                      1            +                                          ω                2                            ⁢                              τ                                  σ                  ⁢                                                                          ⁢                  l                                2                                                                        ∑                  l          =          1                L            ⁢                        ω          ⁢                                          ⁢          τ          ⁢                                          ⁢                      τ                          σ              ⁢                                                          ⁢              l                                                1          +                                    ω              2                        ⁢                          τ                              σ                ⁢                                                                  ⁢                l                            2                                                          Here, τ depends on Q and τσl depends on the specific frequency range. Because the τ method cannot simulate compensation effects, traditional methods using the τ method are inexact.        
Conventional Q-RTM computer methods for imaging seismic data (hereinafter, “Q-RTM computer methods”) require a substantial amount of computer memory because such conventional computer methods use a global approach. The global approach employs a global differential operator to compensate the effects of anelastic attenuation in the conventional Q-RTM computer methods. One global differential operator is a pseudo-spectral operator based on a Fast Fourier Transform. Implementing a Fast Fourier Transform requires the use of information from the entire computational domain. In the global approach, as each seismic wave measurement is recursively traced backward in time from the receiver to the subsurface reflection point, at each backward time step, the calculation of the back-propagated seismic wavefield at each subsurface point requires the wavefield information for the entire computational domain.
The dependency on a global differential operator for each calculation and the size of the global differential operator cause pseudo-spectral computer methods for imaging seismic data (hereinafter, “pseudo-spectral computer methods”) to be computationally expensive, inefficient, and not suitable for parallelization. Computer systems using one or more computing nodes and/or computer clusters may be used for Q-RTM computer methods. Pseudo-spectral computer methods used as part of the process of imaging seismic data often result in poor performance by the computer system on which the conventional Q-RTM computer methods are processed.
Conventionally, for large-scale seismic imaging, due to CPU or graphic processing unit (GPU) memory limitations, the entire computational domain is decomposed into several computational subdomains. Each CPU or GPU calculates one of the computational subdomains. When the entire computational domain is decomposed, individual CPUs or GPUs share the information for the global differential operator that will be applied. The information and in particular the amount of information communicated between CPUs and GPUs slows computations necessary for the seismic imaging. Pseudo-spectral computer methods either lose capability or lose efficiency because computational domain decomposition is difficult with pseudo-spectral computer methods due to the significantly increased communications between CPUs or GPUs.
More memory-efficient Q-RTM computer methods have not been developed. With conventional computer methods, fractional derivatives in dispersion relation compensation cannot be accurately evaluated by using local numerical differential operators, such as finite difference operators.